Locosselli et al. PNAS modeling analysis
# libraries
library(MCMCglmm); library(coda)
library(parallel)
library(phytools); library(brranching); library(ape)
library(dplyr)Dataset
## 'data.frame': 763 obs. of 18 variables:
## $ zone : Factor w/ 3 levels "Sub-tropics",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ ID : int 1 2 3 4 5 6 7 8 9 10 ...
## $ fam : Factor w/ 53 levels "Anacardiaceae",..: 36 30 25 25 25 25 21 27 30 21 ...
## $ taxagroup : Factor w/ 2 levels "Conifer","Eudicot": 1 2 2 2 2 2 2 2 2 2 ...
## $ species : Factor w/ 281 levels "Abies delavayi",..: 5 267 260 260 260 260 26 44 71 73 ...
## $ longevity : int 311 186 221 523 262 219 243 427 308 123 ...
## $ growth : num 1.94 5.52 4 2.62 NA ...
## $ densityCHAVE: num 0.43 0.56 0.601 0.601 0.601 0.601 0.524 0.638 0.447 0.504 ...
## $ human.influ : int 35 39 56 14 18 22 12 15 12 15 ...
## $ tseas : int 1340 1669 1941 1043 1325 1036 1095 685 1095 685 ...
## $ moistwetq : num 1.18 1.38 2.47 1.96 2.19 1.59 1.32 1.36 1.32 1.36 ...
## $ low.w.rad : num 145 145 136 144 143 ...
## $ temp : num 24.5 22.1 27 22.6 27.7 21.6 25.1 26.5 25.1 26.5 ...
## $ moistdq : num 0.18 0.19 0.09 0.4 0.25 0.52 0.44 0.34 0.44 0.34 ...
## $ cloud : int 302 307 596 207 329 179 157 210 157 210 ...
## $ soil_class : Factor w/ 26 levels "Acrisols","Alisols",..: 13 16 9 19 16 16 9 8 9 8 ...
## $ alt.group : Factor w/ 2 levels "High altitude",..: 1 2 2 2 2 2 2 2 2 2 ...
## $ phylo : Factor w/ 281 levels "Abies_delavayi",..: 5 267 260 260 260 260 26 44 71 73 ...Columns information/metadata:
- zone: global climatic zone
- ID: identity of the observation data
- fam: Family name
- taxagroup: taxonomical group, Conifer or Eudicot
- species: Species name
- longevity: recorde maximum longevity
- growth: recorder growth
- densityCHAVE: wood density
- human.influ: Human influence index
- tseas: temperature seasonality
- moistwetq: water-balance soil moisture index of the wettest quarter of the year
- low.w.rad: lowest weekly radiation
- temp: mean annual temperature
- moistdq:moisture index of the driest quarter
- cloud: cloud cover
- soil_class: soil classification
- alt.group: separation of the populations in “High altitude” and “Lowland”
- phylo: Species names used to match the phylogenetic tree
Standardizing numeric variables:
Phylogeny
The phylogenetic tree was constructed using the Phylomatic software (Webb et al. 2008) based on a revised vascular megatree (Gastauer & Meira-Neto 2016). We estimated the ancestral values of longevity and growth at internal nodes using maximum likelihood under the assumption of Brownian motion for trait evolution (Revel 2013). Auxiliary packages used: phytools, brranching and ape.
Inverse phylogenetic matrix to be used in the model
Model 1
Code of the bivariate phylogenetic linear mixed model with the packages MCMCglmm and coda.
The model takes long yours (or days) to run. The final version was stored to be loaded:
Prior specification for the random effects and residuals
prior <- list(R = list(V = diag(2), nu = 1.002),
G = list(G1 = list(V = diag(2), nu = 1.002), # phylogeny
G2 = list(V = diag(2), nu = 1.002), # species
G3 = list(V = diag(2), nu = 1.002))) # soil_classModel
mod <- MCMCglmm(cbind(log(longevity), log(growth)) ~ trait
+ trait:human.influ
+ trait:densityCHAVE
+ trait:moistwetq + trait:temp + trait:moistdq
+ trait:tseas + trait:low.w.rad + trait:cloud - 1,
random = ~ us(trait):phylo + us(trait):species + idh(trait):soil_class,
rcov = ~ us(trait):units,
family = c("gaussian", "gaussian"),
ginverse = list(phylo = inv.phylo$Ainv),
prior = prior, data = data, pr = TRUE,
nitt = 1000000, burnin = 10000, thin = 200)Model summary
##
## Iterations = 10001:999801
## Thinning interval = 200
## Sample size = 4950
##
## DIC: 1512.182
##
## G-structure: ~us(trait):phylo
##
## post.mean l-95% CI u-95% CI eff.samp
## traitlongevity:traitlongevity.phylo 0.4993 0.19931 0.7953 4950
## traitgrowth:traitlongevity.phylo -0.1955 -0.39760 -0.0275 4950
## traitlongevity:traitgrowth.phylo -0.1955 -0.39760 -0.0275 4950
## traitgrowth:traitgrowth.phylo 0.2062 0.06251 0.3730 4950
##
## ~us(trait):species
##
## post.mean l-95% CI u-95% CI eff.samp
## traitlongevity:traitlongevity.species 0.10901 0.04823 0.17329 4950
## traitgrowth:traitlongevity.species -0.02501 -0.06778 0.01552 4950
## traitlongevity:traitgrowth.species -0.02501 -0.06778 0.01552 4950
## traitgrowth:traitgrowth.species 0.09071 0.04905 0.13286 4739
##
## ~idh(trait):soil_class
##
## post.mean l-95% CI u-95% CI eff.samp
## traitlongevity.soil_class 0.06294 0.02436 0.10914 4950
## traitgrowth.soil_class 0.05488 0.02245 0.09276 5420
##
## R-structure: ~us(trait):units
##
## post.mean l-95% CI u-95% CI eff.samp
## traitlongevity:traitlongevity.units 0.18907 0.16606 0.21423 4950
## traitgrowth:traitlongevity.units -0.05821 -0.07277 -0.04211 4950
## traitlongevity:traitgrowth.units -0.05821 -0.07277 -0.04211 4950
## traitgrowth:traitgrowth.units 0.11571 0.09917 0.13324 4950
##
## Location effects: cbind(log(longevity), log(growth)) ~ trait + trait:human.influ + trait:densityCHAVE + trait:moistwetq + trait:temp + trait:moistdq + trait:tseas + trait:low.w.rad + trait:cloud - 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## traitlongevity 5.172109 4.499176 5.821662 4950 < 2e-04 ***
## traitgrowth 1.129257 0.703709 1.586139 5153 < 2e-04 ***
## traitlongevity:human.influ -0.105172 -0.149308 -0.061976 4950 < 2e-04 ***
## traitgrowth:human.influ -0.023794 -0.064216 0.013438 4310 0.220606
## traitlongevity:densityCHAVE -0.046698 -0.121774 0.031116 4950 0.242020
## traitgrowth:densityCHAVE -0.121748 -0.185407 -0.057769 4950 < 2e-04 ***
## traitlongevity:moistwetq 0.037342 -0.024483 0.100625 4950 0.258182
## traitgrowth:moistwetq 0.038007 -0.015371 0.091839 4950 0.174141
## traitlongevity:temp 0.031938 -0.052077 0.114523 4950 0.454545
## traitgrowth:temp 0.107223 0.033680 0.187452 5202 0.005253 **
## traitlongevity:moistdq 0.148818 0.070339 0.226906 4714 0.000404 ***
## traitgrowth:moistdq -0.035666 -0.104497 0.028710 5462 0.294545
## traitlongevity:tseas 0.049887 -0.037768 0.135077 4950 0.267879
## traitgrowth:tseas -0.073050 -0.149220 0.008314 4950 0.074747 .
## traitlongevity:low.w.rad -0.083994 -0.148693 -0.021337 4950 0.012121 *
## traitgrowth:low.w.rad 0.041644 -0.011022 0.094039 4950 0.113939
## traitlongevity:cloud -0.013544 -0.075064 0.054227 4950 0.685657
## traitgrowth:cloud 0.034173 -0.026544 0.088690 5647 0.239192
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Gelman-Rubin criteria of convergence (4 chains)
4 chains of the model previously ran to make the diagnose of the convergence of the chain:
## Potential scale reduction factors:
##
## Point est. Upper C.I.
## traitlongevity 1 1
## traitgrowth 1 1
## traitlongevity:human.influ 1 1
## traitgrowth:human.influ 1 1
## traitlongevity:densityCHAVE 1 1
## traitgrowth:densityCHAVE 1 1
## traitlongevity:moistwetq 1 1
## traitgrowth:moistwetq 1 1
## traitlongevity:temp 1 1
## traitgrowth:temp 1 1
## traitlongevity:moistdq 1 1
## traitgrowth:moistdq 1 1
## traitlongevity:tseas 1 1
## traitgrowth:tseas 1 1
## traitlongevity:low.w.rad 1 1
## traitgrowth:low.w.rad 1 1
## traitlongevity:cloud 1 1
## traitgrowth:cloud 1 1
##
## Multivariate psrf
##
## 1Trace and density plots for fixed effects
Model 2 lowland
Subsetting data:
Prior specification for the random effects and residuals
prior <- list(R = list(V = diag(2), nu = 1.002),
G = list(G1 = list(V = diag(2), nu = 1.002), # phylogeny
G2 = list(V = diag(2), nu = 1.002), # species
G3 = list(V = diag(2), nu = 1.002))) # soil_classModel
mod.low20 <- MCMCglmm(cbind(log(longevity), log(growth)) ~ trait +
+ trait:temp + trait:temp2 + trait:temp3 - 1,
random = ~ us(trait):phylo + us(trait):species + idh(trait):soil_class,
rcov = ~ us(trait):units,
family = c("gaussian", "gaussian"), ginverse = list(phylo = inv.phylo$Ainv),
prior = prior, data = datalow, pr = TRUE,
nitt = 100000, burnin = 10000, thin = 200)Model’s results previoulsy ran:
Model summary
##
## Iterations = 10001:99801
## Thinning interval = 200
## Sample size = 450
##
## DIC: 1015.704
##
## G-structure: ~us(trait):phylo
##
## post.mean l-95% CI u-95% CI eff.samp
## traitlongevity:traitlongevity.phylo 0.8533 0.32932 1.39428 450
## traitgrowth:traitlongevity.phylo -0.2310 -0.51216 0.01856 450
## traitlongevity:traitgrowth.phylo -0.2310 -0.51216 0.01856 450
## traitgrowth:traitgrowth.phylo 0.2103 0.04624 0.39553 450
##
## ~us(trait):species
##
## post.mean l-95% CI u-95% CI eff.samp
## traitlongevity:traitlongevity.species 0.1316 0.04960 0.23736 450
## traitgrowth:traitlongevity.species -0.0235 -0.08112 0.03525 450
## traitlongevity:traitgrowth.species -0.0235 -0.08112 0.03525 450
## traitgrowth:traitgrowth.species 0.1043 0.05968 0.16160 450
##
## ~idh(trait):soil_class
##
## post.mean l-95% CI u-95% CI eff.samp
## traitlongevity.soil_class 0.07081 0.02655 0.14000 450
## traitgrowth.soil_class 0.05778 0.02368 0.09976 450
##
## R-structure: ~us(trait):units
##
## post.mean l-95% CI u-95% CI eff.samp
## traitlongevity:traitlongevity.units 0.20939 0.17288 0.2472 450.0
## traitgrowth:traitlongevity.units -0.06491 -0.09168 -0.0421 486.5
## traitlongevity:traitgrowth.units -0.06491 -0.09168 -0.0421 486.5
## traitgrowth:traitgrowth.units 0.12699 0.10400 0.1507 679.8
##
## Location effects: cbind(log(longevity), log(growth)) ~ trait + +trait:temp + trait:temp2 + trait:temp3 - 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## traitlongevity 5.2951 4.2999 6.1978 450 <0.002 **
## traitgrowth 1.2774 0.6980 1.7618 450 <0.002 **
## traitlongevity:temp -2.4071 -4.5952 -0.2509 450 0.0267 *
## traitgrowth:temp -0.4205 -1.7737 1.4981 450 0.5867
## traitlongevity:temp2 5.3139 1.8497 8.8137 450 <0.002 **
## traitgrowth:temp2 1.1204 -1.7590 3.8290 450 0.4400
## traitlongevity:temp3 -3.1367 -4.7923 -1.4586 450 <0.002 **
## traitgrowth:temp3 -0.5973 -1.9841 0.7990 450 0.3956
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Trace and density plots for fixed effects
Session Information
## R version 3.6.2 (2019-12-12)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Catalina 10.15.7
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
##
## locale:
## [1] pt_BR.UTF-8/pt_BR.UTF-8/pt_BR.UTF-8/C/pt_BR.UTF-8/pt_BR.UTF-8
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] dplyr_1.0.1 brranching_0.5.0 phytools_0.6-99 maps_3.3.0 MCMCglmm_2.29
## [6] ape_5.3 coda_0.19-3 Matrix_1.2-18 knitr_1.23
##
## loaded via a namespace (and not attached):
## [1] jsonlite_1.6.1 foreach_1.4.4 bold_0.9.0
## [4] gtools_3.8.1 shiny_1.3.2 assertthat_0.2.1
## [7] expm_0.999-4 highr_0.8 animation_2.6
## [10] tensorA_0.36.1 taxize_0.9.8 yaml_2.2.0
## [13] pillar_1.4.6 numDeriv_2016.8-1.1 lattice_0.20-38
## [16] glue_1.4.2 quadprog_1.5-7 phangorn_2.5.5
## [19] uuid_0.1-2 digest_0.6.25 promises_1.0.1
## [22] htmltools_0.3.6 httpuv_1.5.1 plyr_1.8.4
## [25] pkgconfig_2.0.3 httpcode_0.2.0 questionr_0.7.0
## [28] bookdown_0.11 purrr_0.3.4 xtable_1.8-4
## [31] corpcor_1.6.9 later_0.8.0 phylocomr_0.1.4
## [34] cubature_2.0.3 tibble_3.0.3 combinat_0.0-8
## [37] generics_0.0.2 ellipsis_0.3.1 cli_2.0.2
## [40] mnormt_1.5-5 magrittr_1.5 crayon_1.3.4
## [43] mime_0.7 evaluate_0.14 fansi_0.4.1
## [46] nlme_3.1-142 MASS_7.3-51.4 xml2_1.3.2
## [49] tools_3.6.2 data.table_1.12.8 lifecycle_0.2.0
## [52] stringr_1.4.0 plotrix_3.7-8 compiler_3.6.2
## [55] rlang_0.4.7 clusterGeneration_1.3.4 grid_3.6.2
## [58] conditionz_0.1.0 iterators_1.0.10 rstudioapi_0.11
## [61] igraph_1.2.4.2 miniUI_0.1.1.1 rmarkdown_1.13
## [64] codetools_0.2-16 reshape_0.8.8 curl_3.3
## [67] reshape2_1.4.3 R6_2.4.1 zoo_1.8-6
## [70] fastmatch_1.1-0 stringi_1.4.3 rmdformats_0.3.3
## [73] crul_0.8.4 Rcpp_1.0.4 vctrs_0.3.4
## [76] tidyselect_1.1.0 scatterplot3d_0.3-41 xfun_0.7References
Webb C.O., D. D. Akerly, S. W. Kembel, Phylocom: software for the analysis of phylogenetic community structure and character evolution. Bioinformatics, 24, 2098-2100. http://dx.doi.org/10.1093/bioinformatics/btn358. PMid:18678590 (2008).
Gastauer,M., J. A. A. Meira-Neto, An enhanced calibration of a recently released megatree for the analysis of phylogenetic diversity. Brazilian Journal of Biology, 76(3), 619-628 (2016).
Revell, L.J. Two new graphical methods for mapping trait evolution on phylogenies. Methods in Ecology and Evolution 4, 754–759 (2013).